#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <stdint.h>
#include <stdbool.h>

#include "lin_al.h"

#define BUFSIZE 8192

/*
 * Dear diary,
 * this is an implementation of Gaussian elimination
 * algorithm for a common matrix.
 */
extern struct matrix *
gauss(struct matrix *A)
{
	double tmpd;
	int i = 0, j = 0;
	register int p = 0, q = 0;

	/* Building a row echelon form of the matrix */
	while (i < A->m && j < A->n) {
		/* Check if we have a zero as a pivot element */
		if (A->data[i][j] == 0.0) {
			/* Look for a non-zero pivot */
			for (p = i + 1; p < A->m && A->data[p][j] == 0.0; ++p);

			if (p == A->m) {
				/* If we didnt find anything */
				if (j+1 != A->n) {
					/* We dont have a pivot
					 * element in this column.
					 * Next! */
					++j;
					continue;
				} else {
					/* We're at the last col. */
					break;
				}
			} else {
				/* Found a pivot element.
				 * Now we need to swap the rows and
				 * go on further */
				swap_rows(A, i, p);
			}
		}

		if (A->data[i][j] != 1.) {
			tmpd = A->data[i][j];
			for (q = j; q < A->n; ++q) {
				A->data[i][q] /= tmpd;
			}
		}

		/*
		 * X - I * a[i][j] --> X
		 */
		for (p = i + 1; p < A->m; ++p) { /* next row */
			tmpd = A->data[p][j];
			for (q = j; q < A->n; ++q) { /* next column */
				A->data[p][q] =
					(A->data[i][q] * tmpd) - A->data[p][q];
			}
		}
		fputs("---\n", stderr);
		fprint_matrix(A, stderr);
		fputs("---\n", stderr);
		++i;
		++j;
	}

	/*
	 * Transforming a row echelon form into
	 * reduced row echelon form (aka row canonical form).
	 */
	--j; --i;
	while ((i >= 0 && i < A->m) && (j >= 0 && j < A->n))
	{
		if (A->data[i][j] != 1.) {
			--j;
			continue;
		}

		for (p = i - 1; p >= 0; --p) {
			tmpd = A->data[p][j];
			for (q = j ; q < A->n ; ++q)
				A->data[p][q] -= A->data[i][q] * tmpd;
		}
		fputs("===\n", stderr);
		fprint_matrix(A, stderr);
		fputs("===\n", stderr);
		--i;
		--j;
	}

	return A;
}

int
main(int argc, char **argv)
{
	struct matrix	*a;
	unsigned int	i = 0U, j = 0U;

	a = my_malloc(sizeof *a);
	fprintf(stdout, "Enter size of augmented matrix (height, length): ");
	scanf("%d%d", &a->m, &a->n);
	fprintf(stdout, "Enter the augmented matrix:\n\n");
	matrix_calloc(a);
	while (i < a->m) {
		while (j < a->n) {
			(void) scanf("%lf", &(a->data[i][j]));
			++j;
		}
		j = 0;
		++i;
	}
	gauss(a);

	matrix_free(a);
	return 0;
}
